1. Contexto
Por cambio climático entendemos tanto el calentamiento global resultado de las emisiones de gases de efecto invernadero, como el incremento en la frecuencia de las catástrofes naturales causadas por las transformaciones en el clima. Este, representa una amenaza existencial y aunque la comunidad científica lo lleva advirtiendo hace años, los gobiernos anteponen políticas con resultados económicamente favorables a corte plazo a determinadas medidas de reducción del cambio climático.
Existen dos problemas destacados por los que no se están aplicando más medidas para combatir al cambio climático:
- Las tecnologías actuales para combatir el cambio climatico tienen por lo general un coste elevada y no conllevan un impacto economicamente favorable directo.
- Es difícil asegurar que las medidas a favor de erradicar el cambio climático van a afectar por igual a todos lo habitantes del planeta.
Con este proyecto se pretende desarrollar una visión intuitiva de la diferencia entre paises en cuanto al uso de energías renovables y convencionales. También como afectan al cambio climático y los efectos negativos este medidos con la diferencia de temperatura.
Los datos tomados se dividen en tres tablas:
Emissiones\(CO_2\)capita: contiene las emisiones de \(CO_2\) por país y año. Para poder comparar las emisiones anuales entre países, estas se miden en toneladas por persona. Estos datos solo reflejan las emisiones resultantes de la producción dentro del territorio nacional, sin tener en cuenta las emisiones por transporte de bienes, i.e. importaciones/exportaciones. Los datos se han contrastado con los datos de Our World in Data y CDIAC.
ConsumoEnergetico: agrupa el consumo de energías primarias dividido según su fuente (e.g. gas natural, carbón, renovables) por país y año. No todas las variables son estrictamente del mismo tipo ya que, por una parte hay consumos de gas, petróleo y carbón en ExaJoules (\(EJ\)), y por otra parte valores de generación en Tera Vatios por hora (\(TWh\)). Como el objetivo es desarrollar una comparativa entre países, esta diferencia no crea ningún problema. Estas cifras, a diferencia de las emisiones de \(CO_2\), son acumulados por país y no por persona. Los datos se han contrastado con el siguiente informe, BP, Statistical Review of World Energy
CambiosTemperatura: contiene las anomalías de temperatura de la superficie terrestre en grados Celsius, i.e. el cambio de temperatura con respecto a una referencia climatológica que se corresponde con el periodo de 1951 a 1980. Este dataset cubre el periodo de 1961 a 2019 con estadísticas mensuales, estacionales y anuales del cambio de temperatura medio y su desviación estándar respecto al periodo de referencia. Los datos están publicados en Climate Nasa distribuidos por la NASA-GISS.
Codigo_Pais: es una tabla adicional con los nombres de todos los países y sus respectivos códigos siguiendo la ISO reconocida a nivel internacional con la combinación de dos o tres dígitos.
La primera tabla se ha obtenido de CDP, Kaggle y se ha contrastado con las fuentes mencionadas para determinar a que hacen referencia exactamente los datos con el fin de no llegar a conclusiones erróneas. La segunda tabla ha sido obtenida directamente de OWID y la tercera tabla de Temperature Change, Kaggle. La tabla de códigos se encuetra en Country Code, Kaggle.
2. Modelo de datos
2.1 Selección de variables
# Carga de datos
co2 = read_csv("Datos/EmisionesC02capita.csv")
energia = read_csv("Datos/ConsumoEnergetico.csv")
temp = read_csv("Datos/CambiosTemperatura.csv")
cod = read_csv("Datos/Codigo_Pais.csv")
Tras la carga de las tablas se han traducido los nombres de las variables.
La tabla de las temperaturas ha requerido de varios cambios para extraer las anomalías de temperatura. En primer lugar se han seleccionado solo las variables útiles para este informe. En segundo lugar, de la variable Month que contiene valores por mes, estación y año, solo se han tomado los valores anuales. La variable Element contiene la desviación estándar y el valor puntual de temperatura por lo que se ha filtrado solo el segundo valor. A continuación, se han pivotado las columnas en las que aparecen los años como variables. Los años se han usado para crear una nueva variable (Año) y los valores de temperatura se han añadido a una columna adicional (Dif_Temp). Esta variable representa la diferencia de temperatura para un año y país concretos con respecto al período de referencia mencionado en le contexto de los datos. Para concluir, se han seleccionadas las variable Pais y Codigo de la tabla de códigos.
Show code
# Selección de variables de cada tabla
co2 %<>% select('Pais' = 'country', 'Codigo' = 'iso_code', 'Año' = 'year',
'CO2' = 'co2_per_capita')
energia %<>% select('Pais' = 'Entity', 'Codigo' = 'Code', 'Año' = 'Year',
'Carbon' = 'Coal Consumption - EJ', 'Gas' = 'Gas Consumption - EJ',
'Petroleo' = 'Oil Consumption - EJ', 'Nuclear' = 'Nuclear Generation - TWh',
'Biomasa' = 'Geo Biomass Other - TWh', 'Hidraulica' = 'Hydro Generation - TWh',
'Solar' = 'Solar Generation - TWh', 'Eolica' = 'Wind Generation -TWh')
temp %<>% select('Pais' = 'Area', Months, Element, Y1961:Y2019) %>%
filter(Months == 'Meteorological year' & Element == 'Temperature change') %>%
pivot_longer(c('Y1961':'Y2019'), names_to = 'Año', values_to = 'Dif_Temp') %>%
select(Pais, Año, Dif_Temp)
temp$Año = as.numeric(gsub('Y', '', temp$Año))
cod %<>% select('Pais' = 'Country_name', 'Codigo' = 'code_3digit')Las tres tablas tienen distintos rangos de años, siendo el registro relacionado con la energía el más nuevo por lo que ha marcado el rango de años de las otras dos tablas.
También es diferente el número de países. Con un rápido análisis se ha observado como la tabla de energías solo contiene información de países con una población relevante, mientras que las otras dos monitorizan valores de todo el mundo. Por ejemplo, no se tiene el consumo de energía en la Antártida pero si su temperatura y \(CO_2\). El elevado número de valores de está variable también se debe a que existen valores que agrupan países y a su vez están listados los países que se encuentran en dichas zonas. Es para unificar el nombre de los países, o áreas territoriales, y evitar agrupaciones que se ha añadido la tabla de códigos.
Show code
# Rango de Años y Países
resumen = data.frame(Tabla = c("Energía", "CO2", "Temperatura"),
"Inicio_Registro" = c(min(energia$Año), min(co2$Año), min(temp$Año)),
"Fin_Registro" = c(max(energia$Año), max(co2$Año), max(temp$Año)),
"Países" = c(length(unique(energia$Pais)), length(unique(co2$Pais)), length(unique(temp$Pais)))
)
kable(resumen, caption = 'Rango de Años y Total de paises', booktabs = T) %>%
kable_styling(latex_options = 'hold_position', font_size = 14)| Tabla | Inicio_Registro | Fin_Registro | Países |
|---|---|---|---|
| Energía | 1965 | 2019 | 100 |
| CO2 | 1750 | 2019 | 236 |
| Temperatura | 1961 | 2019 | 284 |
Para pode analizar la relación entre los tres factores (energía, temperatura y \(CO_2\)) en los distintos países, se han unido las tablas mediante un inner join. De esta manera los datos de los países y años resulten estarán completos. En primer lugar, se ha añadido la variable de códigos a las temperaturas y, a continuación, se han añadido las variables de energía y \(CO_2\) mediante la Key: (Código, Año). En segundo lugar, se han seleccionado las variables de interés y, finalmente, se han ordenado las observación manteniendo primero el orden alfabético del nombre del país y luego por orden cronológico.
# Inner Join
datos <- temp %>% inner_join(cod, by='Pais') %>%
inner_join(energia, by=c('Codigo','Año')) %>%
inner_join(co2, by=c('Codigo','Año')) %>%
select('Pais', 'Codigo', 'Año', 'Dif_Temp','CO2','Carbon':'Eolica')
# Orden cronológico
datos <- datos[order(datos$Año, decreasing=F),] %>% as.data.frame()Show code
# Tabla Interactiva 1
datos %>% mutate(across(where(is.numeric), ~round(.,3))) %>%
mutate(across(where(is.character), as.factor)) %>%
datatable(rownames = F,
options = list(pageLength = 5),
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;',
'Tabla Interactiva 1: ',
htmltools::em('Resultado de la unión de las tablas')
))
Para concluir la selección de variables se ha cambiado el formato de la información de las energías renovables para contar con una variable categórica. Teniendo en cuenta que cada observación es una combinación de país y año, se han creado las siguientes variables:
- Renovable: Fuente de energía renovable que más energía generó. Las observaciones cuya generación total por fuentes renovables es nula se ha indicado con el valor NO.
- Gen_Ren: Total de energía generada a partir de fuentes renovables, en TWh.
A continuación, se han eliminado las variables de generación energética por fuentes renovables.
Show code
# Nuevas variables
renovables = colnames(datos)[10:13]
ren_fun = function(row){
if (sum(row)!=0) {
tipo = which.max(row)
return(renovables[tipo])
}
else {
return("NO")
}
}
# Renovable
datos$Renovable = apply(datos[,10:13], 1, FUN = ren_fun)
# Gen_Ren
datos %<>% mutate(Gen_Ren = Biomasa + Hidraulica + Solar + Eolica,) %>%
select(-c("Biomasa", "Hidraulica", "Solar", "Eolica"))
2.2 Limpieza de datos
Valores Nulos
Observando las 5 primeras filas se detecta que hay dos tipos de valores a tener en cuenta, los que son 0 resultantes del registro de información con un valor de nulo, y los NA, resultantes de no registrar esa información específica.
Show code
| Pais | Codigo | Año | Dif_Temp | CO2 | Carbon | Gas | Petroleo | Nuclear | Renovable | Gen_Ren |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 207 | 0 | 0 | 0 | 6 | 0 | 0 | 0 |
Analizando los países y períodos con NA’s, se ha encontrado una relación entre los períodos de valores faltantes con el resultado de situaciones históricas. Con la disolución en el 1991 de la Unión Soviética y en el 1992 de la República Federal Socialista de Yugoslavia, surgieron varios países como Azerbayan, Estonia o Rusia por un lado y Croacia o Eslovenia de la otra división. Esto explica que haya 12 países con NA’s en la variable de temperatura hasta el 1992. De forma similar Eslovaquia presenta NA’s hasta el 1993 cuando se disolvió Checoslovaquia. Entonces estos hechos históricos tienen un efecto sobre una de las variables, Dif_Temp. Eliminando estas observaciones se perdería la información de las otras variables que si está registrada entonces se han sustituido los NA’s por 0’s, para conservar las otras variables. De ser necesario esto se tendrá en cuenta en futuros análisis.
Otra opción sería analizar los países que hacen frontera con estos problemáticos y hacer una media ponderada de sus valores de temperatura, pero eso requeriría un análisis más exhaustivo teniendo en cuenta la dimensión de cada país y sus respectivas fronteras.
Show code
[1] "Belgium" "Luxembourg" "Slovakia"
[4] "Singapore" "Trinidad and Tobago" "Azerbaijan"
[7] "Belarus" "Estonia" "Kazakhstan"
[10] "Latvia" "Lithuania" "Russia"
[13] "Turkmenistan" "Ukraine" "Uzbekistan"
[16] "Croatia" "Slovenia" También hay NA’s en el registro de petróleo de Bangladesh hasta el 1971, cuando se independizo de Pakistán. Las variables de energía no representan valores percápita, si no el total del país por lo que no sería correcto asignar los valores de Pakistán a estos fuera de registro. Del mismo modo que han procedido al registrar los datos de las variables de energía de Bangladesh para ese mismo período de tiempo, se han asignado a un valor de 0.
A raíz del descubrimiento de valores no registrados se ha procedido a comprobar que todos los países tengan el mismo rango de años.
P65 = datos %>% filter(Año == 1965) %>% select(Pais)
P19 = datos %>% filter(Año == 2019) %>% select(Pais)
setdiff(P19$Pais,P65$Pais) [1] "Azerbaijan" "Belarus" "Croatia" "Estonia"
[5] "Kazakhstan" "Latvia" "Lithuania" "Russia"
[9] "Slovenia" "Turkmenistan" "Ukraine" "Uzbekistan" Las observaciones de estos países no empiezan en el 1965 como el resto.
Show code
P85 = datos %>% filter(Año == 1985) %>% select(Pais)
URSS <- setdiff(P85$Pais, P65$Pais) #1985 aparecen 10 países de la URSS
P90 = datos %>% filter(Año == 1990) %>% select(Pais)
RFSI <- setdiff(P90$Pais, P85$Pais) #1990 aparecen dos países de la Rep.Fed.Soc. de YugoslaviaSe ha confirmado que se empezaron a tener datos de Croacia y Eslovenia en el 1990 y en el 1985 el resto de países.
Valores Atípicos
Se ha considerado que debido a las diferencias de población y territorio entre países, los outliers deben ser detectados por países.
# Detección de Outliers
out_det <- function(data, paises, años, n){
outliers = list()
for ( pais in paises ) {
p = data %>% filter(Pais == pais)
outliers[[pais]] = hampel(p[,n], k = 10, t0 = 3)$ind #Longitud ventana = 2*k+1
}
if (isempty(outliers)) {
return("0 Outliers")
}
else{
df = as.data.frame(plyr::ldply(outliers, rbind))
colnames(df)[1] = "Pais"
df %<>%
pivot_longer(!Pais ,values_to = "indx") %>%
select(Pais,indx)
for (i in 1:nrow(df)){
if (df$Pais[i] %in% URSS){
df$indx[i] = años[df$indx[i + 20]]
}
if (df$Pais[i] %in% RFSI){
df$indx[i] = años[df$indx[i + 25]]
}
else {
df$indx[i] = años[df$indx[i]]
}
}
return(na.omit(df))
}
}
paises = unique(datos$Pais)
años = unique(datos$Año)La función out_det() implementa un hample filter para la detección de valores atípicos en series temporales y hace un loop para analizar cada país por separado. Esta se ha aplicado a cada una de la variables.
Show code
# Se ha procedido igual con cada variable
# cada una con su respectiva n, i.e columna
outliers_temp = out_det(datos,paises, años, n = 4)Todas las variables excepto Gas tienen valores atípicos. A modo de ejemplo se muestran los outliers del consumo de petróleo:
Show code
# Outliers Consumo Petróleo
colnames(outliers_petr) = c("País", "Año")
kable(outliers_petr, caption = 'Outliers Petróleo', booktabs = T) %>%
kable_styling(full_width = F,latex_options = 'hold_position', font_size = 14) %>%
column_spec(1,width = "3cm") %>%
column_spec(2,width = "3cm")
| País | Año |
|---|---|
| Finland | 1995 |
| Kuwait | 1990 |
| Kuwait | 1991 |
| Sweden | 1979 |
Todos los valores atípicos han sido analizados y ninguno de ellos es resultado de una medición errónea de las variables. Entonces, se ha decidido no reemplazarlos ya que representan cambios significativos reales que afectan a los países.
La siguiente gráfica muestra la frecuencia de valores atípicos por año.
Show code
# Frecuencia Outliers
outliers <- rbind(outliers_temp, outliers_co2, outliers_carb,
outliers_petr, outliers_nuc, outliers_renov)
outliers %>% mutate(Año = as.factor(indx)) %>%
ggplot(aes(x = Año, fill = Año, color = "tomato1")) +
geom_bar(stat = 'count', alpha = 0.8) +
My_Theme +
theme(axis.text.x = element_text(angle = 90), legend.position = "none") +
ggtitle("Frecuencia de Outliers", subtitle = "Por Año")
Como futuros proyectos se podría realizar un análisis con más profundidad de los valores atípicos teniendo en cuenta si son anómalos por encima o por debajo del resto de valores y ver si están relacionados los valores atípicos entre variables.
2.3 Tipología de datos
Show code
# Tipología de datos
datos %<>% mutate(across(Año,as.integer), across(c(Pais,Renovable), as.factor))
str(datos)'data.frame': 3710 obs. of 11 variables:
$ Pais : Factor w/ 72 levels "Algeria","Argentina",..: 1 2 3 4 6 8 9 10 11 12 ...
$ Codigo : chr "DZA" "ARG" "AUS" "AUT" ...
$ Año : int 1965 1965 1965 1965 1965 1965 1965 1965 1965 1965 ...
$ Dif_Temp : num -0.135 0.063 0.127 -0.882 -0.196 0 -0.079 -0.461 -0.872 0.119 ...
$ CO2 : num 0.525 2.654 10.683 5.221 0.056 ...
$ Carbon : num 0.00293 0.03044 0.7273 0.21281 0 ...
$ Gas : num 0.02675 0.148442 0.000106 0.063103 0 ...
$ Petroleo : num 0.0555 0.9365 0.652 0.2351 0 ...
$ Nuclear : num 0 0 0 0 0 0 0 0 0.128 0 ...
$ Renovable: Factor w/ 5 levels "Biomasa","Eolica",..: 3 3 3 3 4 3 3 3 3 3 ...
$ Gen_Ren : num 0.4 1.23 7.99 16.08 0 ...- Factor: Las dos variables tipo factor, son categóricas que describen el país y el tipo de tecnología renovable que más energía generó. Los factores facilitan agrupaciones y aplicaciones de modelos estadísticos.
- Character: El código del país recibe un uso como etiqueta por lo que se ha mantenido como char.
- Number: Las otras variables son numéricas siendo el año un integer y las demás de tipo doble ya que pueden tener decimales.
3. Análisis Exploratorio
En este apartado se ha analizado la distribución de las variables agrupadas por países. A continuación, se ha entrado con más detalle a analizar los valores relacionados con el consumo energético y finalmente se ha analizado la relación de las diferéncias de temperatura con las emisiones de \(CO_2\) percápita.
# Estadísticos
stats <- datos %>%
pivot_longer(cols = c(Dif_Temp, CO2, Carbon, Gas, Petroleo, Nuclear, Gen_Ren),
names_to = "Variable", values_to = "Valores") %>%
group_by(Pais, Variable) %>%
summarise(Media = mean(Valores), Desv.Est. = std(Valores),
Min = min(Valores), Q1 = quantile(Valores, 0.25),
Mediana = median(Valores), Q3 = quantile(Valores, 0.75),
Max = max(Valores), Sumatorio = sum(Valores))
Show code
# Tabla interactiva 2
stats %>% mutate(across(where(is.numeric), ~round(.,3))) %>%
mutate(across(where(is.character), as.factor)) %>%
datatable(rownames = F, filter = 'bottom',
options = list(pageLength = 5),
caption = htmltools::tags$caption(
style = 'caption-side: bottom; text-align: center;',
'Tabla Interactiva 2: ',
htmltools::em('Distribución de las variables contínuas por país, datos del 1965 hasta 2019')
))
Por una parte, la media y la desviación típica dan una intuición de como se comporta la función de distribución de las variables. Los valores mínimo, máximo y los cuartiles marcan la distribución del rango de valores de las variables. Y para acabar el sumatorio aporta una medida de la dimensión de la variable durante todo el período del 1965 al 2019. El interés de estos valores en concreto radica en poder comparar los países con la tabla interactiva.
3.1 Consumo energético
Para poder visualizar conjuntamente la energía que proviene de distintas fuentes, se han transforman las energías nuclear y renovables a \(EJ\), donde \(1000TWh = 3.6EJ\).
Show code
# Consumo por país
cons_por_pais19 = datos %>%
filter(Año == 2019) %>%
mutate(Nuclear = Nuclear*3.6*10^-3,
Renovables = Gen_Ren*3.6*10^-3,
`Consumo Total` = Carbon + Gas + Petroleo + Nuclear + Renovables)
# Top 10
top_10 = cons_por_pais19 %>%
arrange(-`Consumo Total`) %>%
head(10) %>%
pivot_longer(c( "Consumo Total", 'Carbon':'Nuclear', "Renovables"),
names_to = 'Fuente',
values_to = 'Consumo') %>%
mutate(Pais = factor(Pais, levels=c(c("China", "United States", "India", "Russia", "Japan",
"Germany", "Canada", "Saudi Arabia", "Brazil", "Indonesia"))))
top_10 %>%
ggplot(aes(x = Pais, y = Consumo, fill = Pais)) +
geom_bar(stat = "identity", alpha = 0.8) +
facet_wrap(~ Fuente, scales = "free") +
ylab("Consumo (EJ)") + xlab("País") +
My_Theme2 +
theme(axis.title.x = element_blank(), axis.text.x = element_blank(),axis.ticks.x = element_blank()) +
ggtitle("Países con mayor consumo energético", subtitle = "Top 10, 2019, por fuente energética primaria")
El conjunto de gráficos muestra como China y EE.UU son los mayores consumidores de energía aunque dependiendo de la fuente primaria hay algunas variaciones en el ranking. El dato más destacable es el caso de Rusia que consume más energía producida con gas natural que China y el hecho que Arabia Saudi se encuentre en el top 10 de consumo total basando su consumo tan solo en petróleo y gas natural.
# Evolución Consumo Total
consumo_total = datos %>%
group_by(Año) %>%
summarise(Carbon = sum(Carbon), Gas = sum(Gas), Petroleo = sum(Petroleo),
Nuclear = sum(Nuclear)*3.6*10^-3, Renovables = sum(Gen_Ren)*3.6*10^-3)Show code
consumo_total %>%
pivot_longer(c('Carbon':'Renovables'),
names_to = 'Fuente',
values_to = 'Consumo') %>%
mutate(Fuente = factor(Fuente, levels=c("Renovables","Nuclear","Gas","Petroleo","Carbon"))) %>%
ggplot(aes(x = Año, y = Consumo, fill = Fuente)) +
geom_area(alpha = 0.8) +
scale_x_continuous(breaks = seq(1965,2020,5)) +
ylab("Consumo(EJ)") +
ggtitle("Consumo Energético Total(EJ), 1965 a 2019") +
My_Theme -> plot3
ggplotly(plot3)En este gráfico se puede ver como el consumo energético total se ha más que cuadriplicado en las últimas 6 décadas, especificamente para los países contenidos en estos datos, de 118.671 a 477.48 Exajoules. Este incremento coincide con el aumento de la población mundial que durante este período, se ha duplicado, y con que la mayoría de países han aumentado sus niveles de consumo.
# Evolución Consumo proporcional
consumo_prop = consumo_total %>%
mutate(Año, sum = Carbon + Gas + Petroleo + Nuclear + Renovables,
Carbon = 100*Carbon/sum,Gas = 100*Gas/sum,
Petroleo = 100*Petroleo/sum,Nuclear = 100*Nuclear/sum,
Renovables = 100*Renovables/sum) %>% select(-c(sum))El siguiente gráfico representa la evolución de la proporción de energía consumida de todos los países, según la fuente primaria.
Show code
consumo_prop %>%
pivot_longer(c('Carbon':'Renovables'), names_to = 'Fuente', values_to = 'Consumo') %>%
mutate(Fuente =factor(Fuente, levels=c("Renovables","Nuclear","Gas","Petroleo","Carbon"))) %>%
ggplot(aes(x = Año, y = Consumo)) +
geom_line(aes(color=Fuente),size = 1.5, alpha=0.7) +
scale_x_continuous(breaks=seq(1965,2020,5)) +
ylab("Consumo%") +
ggtitle("Proporción de consumo energético por fuente primaria",
subtitle = "Entre los años 1965 y 2019")+
My_Theme -> plot4
ggplotly(plot4)En general el uso del petróleo ha disminuido en las últimas 3 décadas. El uso de gas natural como fuente de energía ha aumentado y el del carbón ha disminuido con respecto a los años 60. La energía nuclear y la que proveniene de fuentes renovables siguen siendo las más bajas aunque han aumentado un 6.66% y un 153.1% respectivamente en los últimos 20 años.
3.2 Temperatura y CO2
Varios estudios muestran la relación entre las emisiones de \(CO_2\) y el aumento de la temperatura global de la tierra. Pero, ¿hay relación entre las emisiones de \(CO_2\) percápita de un país y los cambios de temperatura que sufre?
Los siguientes valores son el promedio entre todos los países del conjunto de datos.
Show code
datos %>% group_by(Año) %>%
summarise(meanco2 = round(mean(CO2),2), meantemp = round(mean(Dif_Temp),2)) %>%
ggplot(aes(y=meantemp, x=Año)) +
geom_point(aes(color=meanco2), alpha = 0.8, size=3.5) +
geom_line(color="#000099") +
geom_hline(yintercept = 0, linetype = "dashed", color = "#FF0000") +
annotate(geom = "text", y = -0.2, x = 2010, label = "Respecto el período del 1951 al 1980", color="#FF0000") +
scale_colour_gradientn(colours=rainbow(4)) +
scale_x_continuous(breaks=seq(1965,2020,5)) +
labs(title = "Evolución de la Diferencia de Temperatura promedio",
y = "Temperatura (ºC)",
color = "Emisiones per capita (kg)") +
My_Theme -> emisiones
ggplotly(emisiones) Se puede observar como hasta mediados de los 80 la temperatura era más o menos estable. A partir de los 90, a excepción de algunos años esta temperatura no ha hecho más que aumentar.
# CO2 Vs Diferencias de temperatura
datos %>% group_by(Pais) %>%
summarise(CO2 = round(mean(CO2),3), Dif_Temp = round(mean(Dif_Temp),3)) %>%
ggplot(aes(y=CO2, x=Dif_Temp, color = Pais)) +
geom_point(size= 4, alpha = 0.8) +
theme(legend.position= "none") +
labs(title="CO2 Vs Diferencia Temperatura") +
My_Theme -> co2_vs_temp
ggplotly(co2_vs_temp)
No existe ninguna relación directa entre las anomalías de temperatura y las emisiones de \(CO_2\) por persona. Analizando la distribución de los países respecto al eje horizontal parece que la diferencia media de temperaturas está relacionado con la ubicación geográfica de los países. En el siguiente enlace hay un gráfico que muestra la evolución del cambio de de temperaturas en la superficie terrestre sobre un mapa mundial y se observar como se desplaza en forma de varías franjas de temperatura que cubre varios países simultaneamente.
4. Aprendizaje Estadístico
En este apartado se han aplicado métodos de aprendizaje estadístico a los datos analizados previamente. Inicialmente se ha realizado un análisis de componentes principales que junto con el clustering han representado dos formas distintas de agrupar los países. En segundo lugar, se han comparado el algoritmo SVM y el Random Forest para tratar de clasificar los países según la fuente de energía renovable con la que más energía generan. Finalmente, se ha comparado la eficiencia de una RLM implementando un método selección de variables (Forward Stepwise) con la de un Random Forest, esta vez en la predicción de las variaciones de temperatura de la superficie terrestre.
4.1 Analisis de Componentes Principales
Para el siguiente análisis se han escogido los datos del año más reciente, 2019, para ser más próximo a la realidad actual.
# ACP
datos %>% filter(Año == max(Año)) %>%
select(-c(Pais, Codigo, Año)) -> datos_acp
paises = datos %>% filter(Año == max(Año)) %>% select(Codigo)
rownames(datos_acp) = paises$Codigo
# ONE-HOT Encoding
dummy = dummyVars(~ Renovable, data = datos_acp )
Renovable_dummy = as.data.frame(predict(dummy, newdata = datos_acp))
Renovable_dummy %<>% select(c(Biomasa = Renovable.Biomasa, Eolica = Renovable.Eolica,
Hidraulica = Renovable.Hidraulica, Solar = Renovable.Solar))
# Ensamblaje
datos_acp %<>% cbind(Renovable_dummy) %>% select(-c(Renovable))
cp <- prcomp(datos_acp, scale = T)Show code
#Variación Explicada Acumulada
cp.var = as.data.frame(cp$sdev^2) #Varianzas
cp.expl = cp.var/sum(cp.var) #proporción de variación explicada
colnames(cp.expl) = "var_explicada"
ggplot(data = cumsum(cp.expl), aes(var_explicada, x=1:11))+
geom_line(color='green2') + geom_point() +
geom_hline(yintercept=1, linetype='dashed', color= "#FF3333") +
geom_bar(stat = "identity", alpha = 0.8, fill = "#99FFFF")+
scale_x_continuous(breaks=seq(1,11,1)) +
ggtitle('Variación Explicada Acumulada') +
xlab('Componentes Principales') + ylab('Proporción de Variación') +
My_Theme2
Se puede afirmar entonces que las 4 primeras componentes principales explican más de un \(75\%\) de la variabilidad de los datos. Respecto estas componentes, en la tabla 4 se pueden ver sus cargas que describen la relación entre las respectivas componentes y las variables originales.
- CP1: se aproxima a un indice negativo del consumo energético de los países.
- CP2: representa las emisiones de \(CO_2\) per cápita y el tipo de fuente renovable predominante con una mayor contribución por parte de la hidráulica.
Las otras componentes ya son una mezcla en la se observa parte de la contribución de la temperatura y una alta correlación con tecnologías renovables menos frecuentes que la hidráulica.
Show code
# Cargas
cargas = as.data.frame(cp$rotation[,1:4])
kableExtra::kable(t(cargas), caption = 'Cargas de los 4 primeros componentes principales', booktabs = T) %>%
kable_styling(latex_options = 'hold_position', font_size = 20)
| Dif_Temp | CO2 | Carbon | Gas | Petroleo | Nuclear | Gen_Ren | Biomasa | Eolica | Hidraulica | Solar | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| PC1 | 0.0844983 | -0.0522465 | -0.3882330 | -0.4373038 | -0.4948476 | -0.4403145 | -0.4481719 | 0.0654685 | -0.0580734 | -0.0094363 | 0.0285337 |
| PC2 | -0.0534866 | -0.5014625 | 0.1294046 | -0.1129644 | -0.0487856 | -0.0553438 | 0.1202391 | -0.2283824 | -0.2877482 | 0.6381544 | -0.3925163 |
| PC3 | 0.3615064 | -0.2743393 | -0.1034509 | 0.0273478 | -0.0061710 | 0.1044472 | -0.0618161 | -0.0806431 | 0.6808475 | -0.1783050 | -0.5157356 |
| PC4 | -0.2526690 | -0.0988651 | -0.1608380 | 0.0667643 | 0.0008057 | 0.0250724 | -0.1160086 | -0.8199820 | 0.2513624 | 0.0244250 | 0.3818643 |
En el siguiente gráfico se puede evaluar la estructura de los países. Respecto a la primera componente destacan dos grandes aglomeraciones a la derecha y dos puntos atípicos con valores muy negativos. Tras el análisis de las CP1 y CP2 se entienden los valores de China y EE.UU debido al alto consumo energético de ambas. Luego la separación respecto la segunda componente es menos directa ya que por ejemplo separa Italia de España que producen unas emisiones de \(CO_2/capita\) similares. El hecho por el que se han representado divididas es por la contribución del otro factor de la CP2, es decir Italia emplea la energía hidráulica como su principal fuente de energía renovable, mientras que en España predomina la eólica. Este hecho, explica la mayor parte de la división entre los dos segmentos cerca del eje de ordenadas.
Puede ser interesante realizar consultas con la tabla interactiva 2 para entender mejor esta segmentación de los países.
Show code
# Proyección sobre CP1-CP2
proyec = as.data.frame(cp$x[,1:2])
proyec %<>% setDT(keep.rownames=T)
ggplot(proyec, aes(x = PC1, y = PC2))+
geom_point(size = 1.8, color = "tomato1", alpha = 0.8) +
geom_text_repel(aes(label = rn), size = 2, max.overlaps = 30) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
xlim(-11,4) +
ggtitle("Proyección de los Paises sobre la CP1 y CP2") +
My_Theme2
Show code
#Filtro Códigos
datos %>% filter(Año == max(Año)) %>%
select(Pais, Codigo) %>%
mutate(Codigo = as.factor(Codigo)) %>%
datatable(rownames = F, filter = 'top',
options = list(pageLength = 1, dom = "t"),
caption = htmltools::tags$caption(style = 'caption-side: bottom; text-align: center;',
'Tabla Interactiva 3',
htmltools::em('Filtra el país del cual quieras saber su código o viceversa'))
)Sería interesante para proyectos centrados en ACP, comparar la evolución de las componentes principales, i.e. como cambia el foco de variabilidad entre países con los año. También añadir otras variables para obtener indices que describan el desarrollo de los países o el nivel de adversidad frente a los efecto del cambio climático.
4.2 Clustering de países
K-Means
Siguiendo con la tónica del análisis anterior se ha desarrollado una segmentación de los países usando el algoritmo KMeans.
# KMeans
datos_num = r.datos_acp #datos creados para la acp, con los dummys
from sklearn.cluster import KMeans
sdc = [] #suma de distancias cuadradas al centro del cluster respectivo
for i in range(1,12):
kmeans = KMeans(n_clusters = i, init = "k-means++", max_iter = 300, n_init = 10, random_state = 0)
kmeans.fit(datos_num)
sdc.append(kmeans.inertia_)Show code
# Visualización del error por k
_ = plt.figure()
_ = plt.plot(range(1,12), sdc, "o--")
_ = plt.title("Suma de distáncias cuadradas")
_ = plt.xlabel("Número de Clusters(k)")
_ = plt.ylabel("Suma")
plt.show()
# Ajustar los datos con KMeans
kmeans = KMeans(n_clusters = 4, init="k-means++", max_iter = 300, n_init = 10, random_state = 0)
y_kmeans = kmeans.fit_predict(datos_num)Para aproximarlo a la representación del apartado anterior, se han usado como dimensiones la variable Petroleo(CP1) y la \(CO_2\)(CP2).
Show code
# Visualización de la segmentación
_ = plt.figure()
labels = ["Bajo Consumo de Petróleo", "CHN", "USA", "Consumo Moderado de Petróleo"]
colors = ["red", "blue", "green", "cyan"]
for i in range(0,4):
x = datos_num.Petroleo[y_kmeans == i]
y = datos_num.CO2[y_kmeans == i]
if len(x) > 50:
_ = plt.scatter(datos_num.Petroleo[y_kmeans == i], datos_num.CO2[y_kmeans == i], s = 100, c = colors[i], label = labels[i], alpha=0.2)
else:
_ = plt.scatter(datos_num.Petroleo[y_kmeans == i], datos_num.CO2[y_kmeans == i], s = 100, c = colors[i], label = labels[i])
_ = plt.scatter(kmeans.cluster_centers_[:,4], kmeans.cluster_centers_[:,1], s = 10, c = "yellow", label = "Baricentros")
_ = plt.title("Cluster de países")
_ = plt.xlabel("Consumo Energía fuente Petroleo(EJ)")
_ = plt.ylabel("Emisiones CO2 (kg/capita)")
_ = plt.legend()
plt.show()
Con la regla del codo se ha seleccionado \(k=4\) como el número óptimo de clusters entre los países de estudio. Este valor coincidiendo con la segmentación del ACP. Como en el apartado anterior, se observan China y EE.UU en sus propias burbujas lejanos de cualquier otro país, hecho que no es de extrañar y probablemente aplicable en muchas otras métricas. Con este modelo no se ve tan clara la división de las dos aglomeraciones. A grandes rasgosc el modelo ha vuelto a dividir por la variable Petroleo, y algún otro factor no representado en el gráfico como por ejemplo la variable Hidraulica que con el ACP se ha visto que es un factor importante en la segmentación de los países.
4.3 Clasificación de la fuente de renovable dominante
Los datos de estudio contienen dos variables categóricas, Pais y Renovable. Debido al elevado número de valores de la primera variable, 72 países distintos, se ha seleccionado la segunda con 5 valores como target de la clasificación.
Para entrenar este modelo no se han usado las variables de generación energética de las fuentes renovables. Esto se debe a que la variable dependiente de esta clasificación es función de estas, concretamente el máximo de estas.
Destacar que para este apartado y el siguiente, al tratarse de algoritmos de aprendizaje supervisado se ha usado todo el conjunto de daots, seperando una parte de este para el entrenamiento y otra para el test de los modelo.
Show code
# Valores de la variable objetivo 65-19
datos_clas <- datos %>% select(c(Año, Dif_Temp, CO2, Carbon, Petroleo, Gas, Nuclear, Gen_Ren, Renovable))
frec_ren <- as.data.frame(prop.table(table(datos_clas$Renovable)))
df = data.frame() %>% rbind(round(frec_ren$Freq*100,2))
colnames(df) = frec_ren$Var1
kable(df, caption = 'Distribución de valores(%), 1965-2019') %>%
kable_styling(latex_options = 'hold_position', font_size = 14)| Biomasa | Eolica | Hidraulica | NO | Solar |
|---|---|---|---|---|
| 6.74 | 4.1 | 78.14 | 9.11 | 1.91 |
Para equilibrar un poco las clases se han cogido datos a partir del 2013 que es cuando ya todos los países usan alguna fuente de energía renovable, i.e. desaparece la categoría NO, y las observaciones se distribuyen un poco más en las distintas categorías, como muestra la Tabla 6.
Show code
# Valores de la variable objetivo 13-19
datos_clas %<>% filter(Año>=2013)
frec_ren <- as.data.frame(prop.table(table(datos_clas$Renovable)))
df = data.frame() %>% rbind(round(frec_ren$Freq*100,2))
colnames(df) = frec_ren$Var1
kable(df, caption = 'Distribución de valores(%), 2013-2019') %>%
kable_styling(latex_options = 'hold_position', font_size = 14)| Biomasa | Eolica | Hidraulica | NO | Solar |
|---|---|---|---|---|
| 9.92 | 17.26 | 62.5 | 0 | 10.32 |
Este es un problema de multiclasificación con varias dimensiones. Se han usado dos algoritmos distinos, el SVM, support vector machine y el Random Forest. El siguiente código muestra el preprocesamiento de los datos para ambos modelos.
# Separación de variables
X = r.datos_clas.iloc[:, 0:-1]
y = r.datos_clas.iloc[:,-1]
#Train/Test split (75/25)
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.25, random_state = 0)
#Estandarización
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)Support Vector Machine
Se han entrenado dos modelos de SVM uno con kernel rbf y otro polinomial de grado 3. Seguidamente se ha medido su eficiencia clasificando el set de entrenamiento mediante dos métricas, el porcentaje de observaciones correctamente clasificadas(accuracy) y el valor f1.
# Ajustar el SVM en el Conjunto de Entrenamiento
from sklearn.svm import SVC
rbf = SVC(kernel = 'rbf', random_state = 1, gamma = "scale")
rbf_pred = rbf.fit(X_train, y_train).predict(X_test)
poly = SVC(kernel = "poly", random_state = 1, gamma = "scale")
poly_pred = poly.fit(X_train, y_train).predict(X_test)
# Evaluación de los modelos
from sklearn.metrics import accuracy_score
from sklearn.metrics import f1_score
from sklearn.metrics import confusion_matrix
rbf_accuracy = accuracy_score(y_test, rbf_pred)
poly_accuracy = accuracy_score(y_test, poly_pred)
# rbf_f1 = f1_score(y_test, rbf_pred, average='weighted')
# poly_f1 = f1_score(y_test, poly_pred, average='weighted')Las dos ultimas lineas de código están comentadas, pero se han ejecutado en una chunk escondida debido a un aviso inevitable de Python de la función f1_score.
Show code
# Precisión
print("", 'Accuracy (RBF Kernel): ', "%.2f" % (rbf_accuracy*100), '\n',
'F1 (RBF Kernel): ', "%.2f" % (rbf_f1*100), '\n',
'Accuracy (Poly Kernel): ', "%.2f" % (poly_accuracy*100),'\n',
'F1 (Poly Kernel): ', "%.2f" % (poly_f1*100)) Accuracy (RBF Kernel): 72.22
F1 (RBF Kernel): 63.41
Accuracy (Poly Kernel): 72.22
F1 (Poly Kernel): 62.90Show code
# Matriz de Confusión SVM
cm = confusion_matrix(y_test, poly_pred)
df_cm = pd.DataFrame(cm, index = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]],
columns = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]])
_ = plt.figure(figsize = (10,7))
_ = sn.heatmap(df_cm, annot=True)
_ = plt.title("Matriz de Confusión")
_ = plt.xlabel("Predicción")
_ = plt.ylabel("Real")
plt.show()
Random Forest
Se ha entrenado un modelo y se ha evaluado de la misma manera que el SVM.
# Ajustar un RF a los datos de train
from sklearn.ensemble import RandomForestClassifier
rf_model = RandomForestClassifier(n_estimators = 500,
criterion = "entropy",
random_state = 0)
# Predicción
rf_pred = rf_model.fit(X_train, y_train).predict(X_test)
# Evaluación
rf_accuracy = accuracy_score(y_test,rf_pred)
rf_f1 = f1_score(y_test, rf_pred, average='weighted')Show code
# Precisión
print("", 'Accuracy (RBF Kernel): ', "%.2f" % (rf_accuracy*100),'\n',
'F1 (RBF Kernel): ', "%.2f" % (rf_f1*100)) Accuracy (RBF Kernel): 92.86
F1 (RBF Kernel): 92.95Show code
# Matriz de Confusion RF
cm = confusion_matrix(y_test, rf_pred)
df_cm = pd.DataFrame(cm, index = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]],
columns = [i for i in ["Biomassa","Eólica","Hidráulica","Solar"]])
_ = plt.figure(figsize = (10,7))
_ = sn.heatmap(df_cm, annot=True)
_ = plt.title("Matriz de Confusión Random Forest")
_ = plt.xlabel("Predicción")
_ = plt.ylabel("Real")
plt.show()
En la precisión de los dos modelos vemos que la del Random forest es significativamente mejor. Esto se ve representado en las matrices de confusión, donde vemos que las dos predicen igual la clase Hidráulica pero el segundo modelo aumenta la precisión en las otras tres clases.
4.4 Predicción de las anomalías de Temperatura
Para poder aplicar algoritmos de predicción a los datos de estudio se ha decidido predecir los valores de Dif_Temp apartir de las otras variables. En la realidad esto no tendría mucho sentido ya que las variables son registradas de forma simultanea por lo que no tiene valor predecir la anomalía de temperatura que ocurre en el mismo momento en que se puede medir.
RLM con selección Forwars Stepwise
# Selección de Variables
datos_reg = datos %>% select(-c(Pais, Codigo))
#ONE-HOT Encoding
dummy_reg = dummyVars(~ Renovable, data = datos_reg)
Renovable_dummy = as.data.frame(predict(dummy_reg, newdata = datos_reg))
Renovable_dummy %<>% select(c(Biomasa = Renovable.Biomasa, Eolica = Renovable.Eolica,
Hidraulica = Renovable.Hidraulica, Solar = Renovable.Solar))
#Ensamblaje
datos_reg %<>% cbind(Renovable_dummy) %>% select(-c(Renovable))
set.seed(12)
#Train/Test 75/25%
train_idx = sample(nrow(datos_reg), nrow(datos_reg)*0.75)
datos.train = datos_reg[train_idx, ]
datos.test = datos_reg[-train_idx, ]
num_var <- as.numeric(ncol(datos_acp) - 1) #Vars explicativas
# Función de predicción
predict = function(object, newdata, id, ...){
form = as.formula(object$call[[2]])
mat = model.matrix(form, newdata)
coefi = coef(object,id=id)
xvars = names(coefi)
mat[,xvars]%*%coefi
}
# 10folds CV
k = 10
folds = sample(1:k, nrow(datos.train), replace=TRUE)
cv.errors = matrix(NA,k,num_var, dimnames =list(NULL , paste(1:num_var)))
for(j in 1:k){
best.fit = regsubsets(Dif_Temp~.,data = datos.train[folds!=j,],
nvmax = num_var, method='forward') #Forward SW
for(i in 1:num_var){
pred = predict(best.fit, datos.train[folds==j,],id = i)
cv.errors[j,i]=mean( (datos.train$Dif_Temp[folds==j]-pred)^2)
}
}
rmse.cv = sqrt(apply(cv.errors, 2, mean))
plot(rmse.cv, pch=19, type="b", main="Error", xlab="Número de variables del modelo")
# Conjunto óptimo de variables
reg.sub = regsubsets(Dif_Temp~.,data = datos.train, nvmax = num_var, method='forward')
coef(reg.sub, id = 4)
(Intercept) Año Petroleo Nuclear Eolica
-61.978279733 0.031379016 -0.012744912 0.000667009 0.210211552 Estas son las 4 variables que forman el modelo óptimo encontrado siguiendo el método de selección forward step-wise. Seguidamente se ha calculado una regresión lineal múltiple con estas variables.
# Modelo Final
reg.f = lm(Dif_Temp ~ Año + Petroleo + Nuclear + Eolica, data = datos.train)
lm.pred = predict(reg.f, newdata = datos.test)
error.rlm = mean((datos.test[, "Dif_Temp"] - lm.pred)^2)
summary(reg.f)
Call:
lm(formula = Dif_Temp ~ Año + Petroleo + Nuclear + Eolica, data = datos.train)
Residuals:
Min 1Q Median 3Q Max
-2.02171 -0.34399 -0.01749 0.32684 1.91119
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.198e+01 1.312e+00 -47.230 < 2e-16 ***
Año 3.138e-02 6.588e-04 47.633 < 2e-16 ***
Petroleo -1.274e-02 3.341e-03 -3.815 0.000139 ***
Nuclear 6.670e-04 1.748e-04 3.815 0.000139 ***
Eolica 2.102e-01 5.354e-02 3.926 8.84e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5285 on 2777 degrees of freedom
Multiple R-squared: 0.4851, Adjusted R-squared: 0.4843
F-statistic: 654 on 4 and 2777 DF, p-value: < 2.2e-16De la regresión anterior todas las variables son estadisticamente significativas al \(1\%\). Respecto la precisión del modelo para explicar las diferencia de temperaturas, el \(R^2_{ajustado}\) indica que la regresión explica alrededor de la mitad de variable dependiente.
Random Forest Regression
Ya que este algoritmo ha sido tan eficiente en la clasificación del apartado anterior, ha sido seleccionado también en la predicción de las diferencias de temperatura. Para entrenar el Random Forest se van han usado las mismas variables y conjuntos de entrenamiento y test que con la RLM.
Mediante un grid searching se ha calculado el valor optimo para el parámetro m del modelo. Este parámetro representa el número de variables, elegidas aleatoriamente, usadas como candidatos en cada split).
Show code
# Tuning de m
mtry = c(2,4,6,8,10)
errors = matrix(NA , length(mtry) , 1, dimnames =list(paste(mtry), NULL))
set.seed(12)
for(i in 1:length(mtry)){
rf = randomForest(Dif_Temp~ . , data = datos.train, mtry = mtry[i],
metric = 'RMSE')
errors[i,1] = mean(rf$rsq)
}
plot(x = rownames(errors), errors, pch=19, type = "b",
main = "Error", xlab= "m", ylab= "RMSE")
Show code
# Modelo Final
rf.f = randomForest(Dif_Temp ~ . ,
data = datos.train,
ntree = 500, mtry = 2,
metric = 'RMSE', importance = T)
# rf.pred = predict(rf.f, newdata = datos.test)
# error.rf = mean((datos.test[, "Dif_Temp"] - rf.pred)^2)
#Importancias de las variables
rf.f$importance %>% as.data.frame() %>%
select(Prop_ECM = "%IncMSE") %>%
cbind(Var = rownames(rf.f$importance)) %>%
ggplot(aes(x = Var, y = Prop_ECM)) +
geom_bar(stat = "identity", aes(fill = Var)) +
theme_minimal() +
theme(axis.title.x = element_blank(), axis.text.x = element_blank(),axis.ticks.x = element_blank()) +
ggtitle('Proporción de reducción Media del ECM')
Un estudio interesante sobre la predicción de anomalías de temperatura que se podría hacer con estos datos es una predicción por serie temporal. En este caso si se estarían prediciendo valores que aún serían desconocidos.
Conclusiones
El objetivo del proyecto era desarrollar una comparativa entre el consumo energético de diferentes países y sus emisiones de \(CO_2\), además de tener una visión de las anomalías de temperatura que estamos sufriendo en los últimos años a nivel global. Creo que con el anterior análisis se ha alcanzado este objetivos. Además, han aparecido nuevas idea más concretas sobre que caminos de investigación pueden ser interesantes en este ámbito. Uno que me interesa en especial es el desarrollo de indicadores para comparar los países. Otro ejemplo sería el análisis de la relación entre las emisiones de \(CO_2\) percápita y la fuente de energía primaria más usada en el país. En conclusión, a mi parecer este ha sido un proyecto muy interesante con el que he aprendido mucho no solo en conceptos, también en intereses profesionales.